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Relaxation dynamics of functionalized colloids on attractive substrates 
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Particle-based simulations are performed to study the post-relaxation dynamics of functionalized 
(patchy) colloids adsorbed on an attractive substrate. Kinetically arrested structures that depend 
on the number of adsorbed particles and the strength of the particle-particle and particle-substrate 
interactions are identified. The radial distribution function is characterized by a sequence of peaks, 
with relative intensities that depend on the number of adsorbed particles. The first-layer coverage 
is a non-monotonic function of the number of particles, with an optimal value around one layer 
of adsorbed particles. The initial relaxation towards these structures is characterized by a fast 
(exponential) and a slow (power-law) dynamics. The fast relaxation timescale is a linearly increasing 
function of the number of adsorbed particles in the submonolayer regime, but it saturates for more 
than one adsorbed layer. The slow dynamics exhibits two characteristic exponents, depending on 
the surface coverage. 


I. INTRODUCTION 

Large-scale production of materials with enhanced 
physical properties from the self-organization of colloidal 
particles is believed to set the stage for a revolution in 
materials engineering [iHia. There has been a sustained 
effort towards finding design rules to obtain the desired 
structures through the control of the experimental condi¬ 
tions and the particle-particle interactions. One promis¬ 
ing route is to promote the aggregation of anisotropic col¬ 
loidal particles on an attractive substrate HMH]. How¬ 
ever, the feasibility of the target structures is deeply 
compromised by the kinetic pathways of aggregation and 
post-relaxation. In this work, we study the kinetics of 
relaxation of patchy particles on a flat substrate. 

Patchy particles usually refer to spherical colloids with 
a chemical decoration (patches) on their surface. The 
idea is to control the directionality and strength of the 
particle-particle interaction, the maximum valence of the 
colloid, and the local structure of the aggregates. The 
functionalization of the patches can be obtained using, 
for example, DNA, polymers, or enzymes p!6l[2Q] . This 
type of functionalization allows a very detailed control 
on of colloidal valence, strength, and selectivity of the 
colloidal bonds. During self-assembly, the chemical na¬ 
ture (or strength) of the patch-patch interaction leads 
to large energy barriers impeding the formation of equi¬ 
librium structures, promoting kinetically trapped struc¬ 
tures, such as gels, glasses or polycrystals [DisiHas]. One 
route to control the assembly of non-equilibrium struc¬ 
tures is the use of flat or patterned substrates [26ll29] . 
Recent studies of self-assembly of patchy particles on sub¬ 
strates have focused on the thermodynamic structures 
[3QH34] or on the fully irreversible regime [35H38] . Here, 
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we study the relaxation dynamics on the substrate and 
how it depends on the number of adsorbed particles and 
the strength of the particle-substrate interaction. This 
article is organized in the following way. In the next 
section, the model and the simulation parameters are de¬ 
scribed. Results are discussed in Sec. M and in Sec. IlYl 
we draw some conclusions. 


II. MODEL AND NUMERICAL SIMULATIONS 

We consider a three dimensional system of spherical 
colloidal particles with three patches distributed along 
the equator, with an opening angle of 27r/3. Colloidal 
particles are represented by a core of mass rric and three 
patches at a distance R from the center. The particles 
representing the patches are of zero diameter and mass 
Their relative position to the center of the core 
is fixed at all times. 

The core-core interaction is repulsive, described by a 
Yukawa potential (see Fig. Ba)), 

Vy(r) = Iexp{-k[r-{Ri+Rj)]), (1) 

where Ri and Rj are the effective radii of the two inter¬ 
acting particles, A = 1 is the interaction strength and 
k = 4 the inverse of the screening length. The core-core 
interaction is truncated at a cutoff distance Tc = 3R (at 
Tc the potential is 10~^A/k). The patch-patch interac¬ 
tion is described by an attractive Gaussian potential [25] , 
(see Fig.^a)), 

Voirp) =-€exp{-arl), ( 2 ) 

where e = 40 (in units of ksT) is the interaction 
strength, a = O.IR the width of the Gaussian, and 
is the distance between the two interacting patches. The 
patch-patch interaction is truncated at a cutoff distance 
rpc = 2R. The superposition of these two interactions 
and the resulting energy landscape are represented in 
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FIG. 1. (a) Schematic representation of the core-core re¬ 

pulsive interaction (Yukawa) and the patch-patch attraction 
(Gaussian) as a function of the distance between the center 
of two patchy particles when their patches are aligned, (b) 
Energy landscape of a single-patch particle interacting with 
one single patch as a probe at a position (x,y) relative to the 
center of the patchy particle. 


Figs, [^a) and (b). For these parameters, we expect at 
most one bond per patch. However, as we discuss be¬ 
low, in certain cases we observe the formation of double 
bonds, i.e. two bonds per patch. 

We consider an attractive substrate interacting 
isotropically with the patchy particles. The pairwise po¬ 
tential was derived from the Hamaker theory for two 
spheres [39] in the limit where the radius of one parti¬ 
cle diverges. This gives. 


Ah \ 2R{R + D) ( D \ 

^ 6 [d{D + 2R) \D + 2Rj\' 

for the attractive part and 
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D SR \ 
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( 3 ) 
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for the repulsive one, where Ah is the Hamaker’s con¬ 
stant and D = r — R the distance between the surface 
of the particle and the substrate. The potential is highly 
repulsive for distances shorter than the particle radius 
and is attractive at longer distances. It has a minimum 
at r ^ I.IR. 

To resolve the stochastic trajectories of the particles, 
we describe them as rigid bodies and integrate the cor¬ 
responding Langevin equations of motion for the trans¬ 
lational degrees of freedom. 


mv{t) = -VfV{r) - v{t) + J -^(t), (5) 

n V n 

and rotational ones, 

IQ{t) = -V,V(0) - -iSit) + (6) 

Trp y 7y 

The equations are integrated using the velocity Ver- 
let scheme with a time step of At = 0.01 and Large- 
scale Atomic/Molecular Massively Parallel Simulator 
(LAMMPS) for efficient parallel simulations [40]. v and 
uj are the translational and angular velocity, m and / are 
mass and inertia of the patchy particle, V is the pair¬ 
wise potential, and ^(t) is the stochastic term, from the 
thermal noise, given from a random distribution of zero 
mean. We consider the damping time for the transla¬ 
tional motion. 


m 

dirrjR^ 


( 7 ) 


which we set to be = 0.1 At. From the Stokes-Einstein- 
Debye relation im, 


-Dr 3 

A ^ 4^’ 


(8) 


and so the rotational damping time = lOrt/3. 

The parameters for the Langevin dynamics simula¬ 
tions were taken from Ref. [42]: m = and 

R = 0.5/im. To access time scales of the order of the 
second we used translational and rotational diffusion co¬ 
efficients Dt ~ 6/im^/s and Dr ~ 18<s“^, which are one 
order of magnitude larger than the ones observed exper¬ 
imentally in solution at room temperature [42|. There 
are, at least, two relevant time scales governing the re¬ 
laxation process: The typical time for a bond to break 
and the diffusion time. During our simulations we did 
not observe any bond breaking event, as ksTje is low. 
By contrast, the relaxation dynamics consists of particle 
rearrangements leading to the formation of new bonds. 
Thus, in this low temperature limit, considering a larger 
diffusion coefficient corresponds to re-scaling the time. 

The initial structures were generated using the stochas¬ 
tic model first presented in Ref. [43] where bonds are con¬ 
sidered irreversible during the growth and the mechanism 
of mass transport to the substrate is advective [44]. We 
measure the total number of particles N in monolayers 
(ML) corresponding to particles, where L is the lat¬ 
eral size of the substrate in units of the particle diameter 
2R. 
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N=0.25 


N=1 



FIG. 2. Snapshots for = 0.25 and = 1. Images generated att = {0,10 ^,10 ^,10°,10^} seconds on a substrate of lateral 
size L=32. Colors indicate the number of bonds of the colloidal particle. 


III. RESULTS 
A. Particle bonds 

We analyze the dynamics for N = 0.25 and N = 1. 
From the snapshots in Fig. one can see that for the 
lower N the overall structure evolves from an initially 
homogeneous distribution of patchy particles to a hand¬ 
ful of aggregates. We define aggregates as particles con¬ 
nected through their patches. For N = 1 the spatial 
distribution of particles does not change significantly in 
time and only a single aggregate spanning the entire sub¬ 
strate is formed. This result hints at two different mech¬ 
anisms of relaxation. While for lower N, the particles 
need to diffuse and rotate to maximize the coordination 
of each particle, and thus decrease the overall energy, 
for larger N, the diffusion of the particles is hindered by 
the particle-particle repulsive interaction and relaxation 
is mainly driven by the local rearrangements of individual 
particles. 

In Fig.j^we plot the ratio of particles of a certain num¬ 
ber of bonds as a function of time. A qualitative different 
behavior is observed for the two values of N. The inset 
of Fig. I^a) shows the binding rate, given by the time 
derivative of the total number of bonds. The initial bind¬ 
ing rate for N = 1 is one order of magnitude larger than 
for N = 0.25, as the initial interparticle distance on the 
substrate is lower and particles can promptly form bonds 
with their neighbors. This is also visible from the evolu¬ 
tion of the fraction of isolated particles (unbonded) in the 
main plots of Figs. a) and (b). For N = 1 (Fig. |^b)), 
this fraction vanishes for t ~ 0.1, while for N = 0.25 
Fig. I^a)) it takes six times longer for every isolated par¬ 
ticle to form at least one bond. The fraction of particles 
with one and two bonds is characterized by a maximum 
at an intermediate time. As the initial relaxation dy¬ 


namics is faster for A/" = 1 than N = 0.25, these maxima 
occur much earlier. For longer times, the fraction of par¬ 
ticles with three bonds dominates (at low temperatures 
and strong patch-patch interaction). The fraction of par¬ 
ticles with three bonds is higher for A/" = 1, as only one 
aggregate is formed. Note that, for N=l, about 10% of 
particles have a final number of bonds larger than three 
(the number of patches), while there are almost none for 
N = 0.25. This result suggests that for a sufficiently 
large number of particles on the substrate the interpar¬ 
ticle repulsion is such that multiple bonds per patch are 
formed to decrease the overall energy. 

One expects that three-equally-spaced-patch particles 
self-organize into a honeycomb-like structure [45]. How¬ 
ever, the radial distribution function, shown in Fig. [^a), 
consists of a sequence of peaks at positions that differ 
from the characteristic interparticle distances of the hon¬ 
eycomb lattice. The first peak of the radial distribution 
function is given by particles in direct contact. The sec¬ 
ond one indicates a square-like arrangement of connected 
particles, however they are not periodically reproduced, 
since the third and fourth neighbors of the square lattice 
are not present on the radial distribution function. The 
following peaks of the radial distribution function range 
from a five-particle loop to more stretched ones of six, 
seven, or eight particles. 

Figure [^b) shows the distribution of angles between 
two neighbors connected to one particle (see the inset 
scheme). The first peak indicates the occurrence of dou¬ 
ble bonds on a single patch (see inset snapshot). This 
peak occurs at a = 7r/3 suggesting that the local struc¬ 
ture resembles a triangle. This is an uncommon event 
since the probability that one particle has, at least, one 
double bond is 10% (see discussion above). This event 
has a large effect on the angular distribution as it con¬ 
tributes three times to the statistics, one for each of the 
three particles. These structures are not visible in the ra- 
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FIG. 3. Fraction of particles with a number of bonds i, for 
i — {0,1, 2, 3, > 3}, as a function of time for (a) N — 0.25 and 
(b) iV = 1 on a substrate of lateral size L=64, averaged over 
10 samples. Inset: The binding rate, given by the derivative 
of the total number of bonds (x 10^), as a function of time for 
N — 0.25 and A/" = 1 on a substrate of size L=64. 


dial distribution function of Fig. |^a) since the distance 
between the particles is one diameter. The next peak 
on the angular distribution function is related to square 
structures (7r/2), and the widest peak is the superposition 
of the contribution of the other structures forming loops 
of five to eight particles. This distribution is considerably 
different from the original structure from the stochastic 
model (without relaxation), which exhibits only one peak 
at (a = 27r/3 [43] . 


B. Aggregation dynamics 

We proceed to analyze the aggregation dynamics. We 
label particles belonging to the same aggregate using the 
Hoshen-Kopelman algorithm [46]. As shown in Fig. [^for 
N = 1 the dynamics evolves to form a single aggregate. 
For N = 0.25 the formation of a single aggregate may 
occur at long enough times, if the patch-patch interaction 
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FIG. 4. (a) Radial distribution function A(r, (5r), where r is 
the particle distance and 6r = 0.001 is the size of the bins, for 
N — {0.25,0.5,1} on a substrate of lateral size L=64, aver¬ 
aged over 10 samples, (b) Distribution function of the angle 
a between three connected particles for N = {0.25, 0.5,1} on 
a substrate of lateral size L=64, averaged over 10 samples. 


is strong enough. 

We measure the number of aggregates Ng as a func¬ 
tion of time for different A, as shown in Figs. and 
For a strong patch-patch interaction and very long times 
(and below the bond breaking time), all particles are 
connected, forming a single aggregate. From these re¬ 
sults, we can distinguish two regimes: the first regime 
for t < 0.25 and the second for t >= 0.25. In Fig. 
we plot the number of aggregates as a function of time 
for the first regime to show that Ng has an exponential 
decay, 

Vsi(i) ~ exp(-^it), (9) 

where is the inverse of a characteristic time scale. As 
we can see from the snapshots in Fig. in the first 
regime, the dynamics is dominated by fast formation of 
bonds with local neighbors leading to the formation of 
small aggregates. 

In the inset of Fig. we plot as a function of N. For 
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N=0.25 


N-1 


FIG. 5. Snapshots for N = 0.25 and N = 1. Images generated att = {0,10 ^,10 ^,10°,10^} seconds on a substrate of lateral 
size L=32. Particles of the same color belong to the same aggregate of connected particles. 








FIG. 6. Exponential decay, at initial times, for the number 
of aggregates as a function of time for N from N = 0.2 to 
N = 0.9, on a substrate of lateral size L=64 and averaged 
over 10 samples. Inset: Exponential decay as a function of 

N. 


low N, increases linearly with N. In the submonolayer 
regime, particles are initially distributed on the substrate 
in a homogeneous fashion. They diffuse and form aggre¬ 
gates. The typical area that each isolated particle needs 
to diffuse and aggregate scales with the inverse of the 
total number of particles on the substrate (N). Since 
is the inverse of the characteristic time it should scale 
linearly with N. For N = 0.8 (and above) as increasing 
the number of particles does not reduce the typical inter¬ 
particle distance in the first-layer on the substrate due 
to the nucleation of the second layer, saturates and 
remains practically independent of N. 

Figure. shows the approach to the asymptotic value 
Ns{oo) = 1 for different values of N. We can see a power- 



FIG. 7. Number of aggregates as a function of time for N = 
{0.1,0.2,0.5,0.6,0.7}, on a substrate of lateral size L=64 and 
averaged over 10 samples. 


law behavior at long times, 

7V,2(^)-iV,(oo)~^-«^ (10) 

which reveals the second relaxation regime. As can be 
seen in the snapshots of Fig.[^ the overall structure does 
not change significantly in this regime. From Fig. [^we 
observe two different power laws for small and large N. 
For = 0.1 and N = 0.2, an exponent of ^2 = 0-8 A 
0.1 is found, indicating very slow aggregation dynamics. 
For N > 0.4, we found an exponent of ^2 = 3.4 ± 0.4, 
indicating a faster relaxation. For values of = 0.3 and 
N = 0.4, crossover effects are observed suggesting that 
the transition between the two regimes occurs around 
these values. 

The power-law scaling has one of two characteristic ex¬ 
ponents depending on N. For small N, aggregates need 
to diffuse and eventually merge. While for larger N, dif- 
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FIG. 8. Wrapping probability 11 as a function of N. Results 
are averages over 20 samples for substrates of lateral size of 
L = {16,32,64}. 


N=1 
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FIG. 9. Fraction of particles in the spanning aggregate (Poo) 
at the percolation threshold as a function of the substrate size 
for N = 0.5 (open) and iV = 1 (hlled). Results are averages 
over 10 samples of substrate size ranging from L=4 to L=64. 


fusion is negligible and aggregates need to readjust to 
merge with neighboring ones. This transition can be re¬ 
lated to a percolation transition. In Fig. we plot the 
wrapping probability, given by the probability of finding 
an aggregate that spans the system from one side to the 
other (the same aggregate crosses opposite boundaries 
along the x-direction) as a function of N^ for three dif¬ 
ferent sizes. The percolation transition only occurs for 
A/" > 0.3 which is in agreement with the transition found 
in the second relaxation regime shown in Fig. [^but be¬ 
low the mean-field solution of 0.5 for a Bethe lattice with 
three neighbors m- 

The percolation threshold was estimated for N = 0.5 



FIG. 10. Substrate first-layer coverage as a function of 1/t 
for N — 1 and interaction strengths of the substrate potential 
of Ah = {20,40, 60, 80,100} for a system of linear size L=32 
and averaged over 10 samples. Inset: Asymptotic hrst-layer 
coverage as a function of the substrate potential interaction 
strength for N — 1. 


and = 1 by the peak on the order parameter standard 
deviation. From the position of the peak we can estimate 
the percolation threshold in the thermodynamic limit by 
extrapolating for an infinite substrate. To characterize 
the spanning aggregate we plot in Fig. its size at the 
percolation threshold as a function of the substrate lat¬ 
eral size L for A/" = 0.5 and N = 1. It scales as, 

( 11 ) 

where df is the fractal dimension. For N = 0.5 the fractal 
dimension is df = 1.94 ± 0.02, in agreement with that 
for two-dimensional random percolation. The numerical 
value of the fractal dimension for N = 1 is slightly larger 
than the one expected for two-dimensional percolation, 
due to the contribution of particles from the second layer 
that also belong to the spanning aggregate. 


C. First-layer coverage 

Figure pT] shows the first-layer coverage for A^ = 1 as 
a function of the inverse time. The first-layer coverage is 
defined as. 


e{t) = 


^suhs^col 
Lnr-Lij 


( 12 ) 


where Ngubs is the number of particles on the substrate 
within a distance of one diameter from it, Acoi the cross- 
section area of a sphere, and L the lateral size of the 
substrate. 0 is initially lower than the expected value 
for single-layer random sequential adsorption [48], due 
to the formation of more than one layer m- During 
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FIG. 11. Asymptotic first-layer coverage as a function of N 
for a particle-substrate interaction strength of Ah = 40. In¬ 
set: Magnified region from iV = 1 to iV = 4. Top snapshot: 
Structure formed above de substrate layer that can hinders 
the access of more particles to the substrate. 


the relaxation, 0 increases, suggesting that particles ad¬ 
sorbed on top of others are attracted to the substrate. 
The coverage exceeds the one expected for a Honeycomb- 
like arrangement. We extrapolate the first-layer coverage 
in the thermodynamic limit for different values of the 
strength of the particle-substrate interaction. Ah (see 
inset of Fig. [Tq| ). We observe a dependence on Ah^ as 
we increase the substrate interaction, the local restruc¬ 
turing increases the possibility of more particles being 
adsorbed on the substrate, increasing the first-layer cov¬ 
erage. The radial distribution function for various Ah in¬ 
dicates some inter-particle penetration, a larger peak for 
the square structures, and a more enhanced peak related 
to the contribution of the loops of five to eight particles 
(not shown here). 

In Fig. the extrapolated first-layer coverage as a 
function of N is plotted. We observe the predictable in¬ 


crease in the submonolayer case, since all particles even¬ 
tually move to the first layer. However, for larger number 
of particles, the first-layer coverage decreases with N (see 
inset of Fig. 11). This is related to the formation of con¬ 
nected structures of patchy particles that hold the par¬ 
ticles away from the minimum of the particle-substrate 
potential (see snapshot in Fig. 11). 

IV. CONCLUSIONS 


Self-assembly of strongly interacting particles towards 
thermodynamic structures is typically hindered by the 
formation of kinetically trapped structures that are sta¬ 
ble over long time scales. For the self-assembly of patchy 
colloids on an attractive substrate we considered the re¬ 
laxation dynamics towards equilibrium, at low tempera¬ 
ture, and identified the relevant kinetic structures. We 
have shown how the formation and relaxation of kinetic 
structures depends on the total number of particles and 
on the particle-substrate interaction. We observed a de¬ 
pendence of the relaxation dynamics on the number of 
absorbed particles at both short and long times. We have 
found that the long time dynamics exhibits two charac¬ 
teristic power-law exponents depending on the coverage. 

Combining the directionality of particle-particle inter¬ 
actions and the presence of a substrate provides a promis¬ 
ing route to grow (mono- and multi-layer) regular struc¬ 
tures. However, our results reveal that, while strong 
patch-patch interactions favor the resilience to thermal 
fluctuations they yield large barriers to the relaxation 
towards equilibrium. Thus, as a follow up it would be 
of interest to explore possible annealing strategies (such 
as temperature annealing cycles) to overcome these bar¬ 
riers. Nevertheless, it is interesting to notice that these 
non-equilibrium structures are stable over long periods 
of time. A systematic analysis of their mechanical and 
optical properties might reveal them useful for certain 
practical applications. 
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